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ABSTRACT 

In this paper we report on the improvements implemented in the cosmological radiative 
transfer code CRASH. In particular we present a new multifrequency algorithm for 
spectra sampling which makes use of colored photon packets: we discuss the need for 
the multi-frequency approach, describe its implementation and present the improved 
CRASH performance in reproducing the effects of ionizing radiation with an arbitrary 
spectrum. We further discuss minor changes in the code implementation which allow 
for more efficient performance and an increased precision. 
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1 INTRODUCTION 

The proper and accurate treatment of radiative transfer 
(RT) in cosmological numerical simulations is a long stand- 
ing problem, which has received great and ever increasing 
attention in the last decades. Despite this, the progress done 
in numerical RT has been somewhat slower with respect to 
other fields of numerical cosmology and astrophysics, like 
Large Scale Structure simulations or Magneto Hydrodynam- 
ics schemes, which in the last years underwent an impressive 
series of successes and refinements. The reason for the lag 
above is mainly due to the extreme complexity underlying 
the process of radiative transfer, which in a cosmological 
setting couples also with gas dynamics, cosmic expansion, 
structure formation and chemistry. 

The seven dimension cosmological RT equation does not 
allow an exact solution, neither analytic nor numerical, and 
the need to deal with arbitrary geometries and a huge dy- 
namical range in the optical depth distribution prevents the 
adoption of simplifying assumptions to reduce the dimen- 
sionality of the equation and the complexity of the problem. 
Several methods have been proposed in the literature, re- 
lying on a variety of assumptions and approximations, and 
often optimizing the treatment of RT for different problems, 
namely different astrophysical configurations. For a detailed 
review of the various methods developed in the literature we 
refer the reader to Iliev et a/. (2006; 106), where a system- 
atic comparison of 11 cosmological radiative transfer codes 
through a set of tests is described. Since the publication of 
the paper, more RT codes designed for cosmological applica- 
tions have appeared in the literature (Altay, Croft & Pelu- 
pessy 2008; Pawlik & Schaye 2008; Finlator, Ozel & Dave 
2008). 

The strategy adopted by our group to face the problem 
of cosmological radiative transfer is based on a ray-tracing 



Monte Carlo (MC) scheme which applies to 3D Cartesian 
grids and which is based on the particle nature of the ra- 
diation field. This approach offers several advantages. The 
radiation is described in terms of photons which, grouped 
for computational convenience into monochromatic photon 
packets, travel through the simulation volume along rays. 
This allows to solve the transfer in one dimension: (i) hy 
following monochromatic packets of photons along rays one 
gets rid of the explicit dependence on direction and fre- 
quency, and (a) by casting rays into the assigned grid, one 
can furthermore remove the dependency on the position. 
Instead of solving directly for the intensity of the radia- 
tion field, /i/(f, f2,t), one just needs then to model the in- 
teraction of the photons with the gas inside a cell. This 
approach allows to enormously simplify the treatment of 
radiative transfer, and to leave the geometry of the prob- 
lem completely arbitrary. On the other hand, as drawback, 
ray tracing techniques are usually computationally expen- 
sive and tend to suffer from low angular resolution. These 
problems are partially overcome by adopting a statistical 
description of the radiation field, by Monte Carlo sampling 
directions and spectral energy distributions. Differently from 
short and long characteristic methods, in which the rays are 
cast along fixed directions, in MC ray tracing the rays are 
shot in random directions, allowing for a more efficient im- 
plementation of multiple sources and a diffuse component of 
the ionizing radiation, e.g. recombination radiation or the ul- 
traviolet background (UVB), as well as anisotropic angular 
emissivities. These advantages come with the introduction 
of numerical noise, intrinsic to all MC schemes, which can 
be nevertheless kept arbitrarily low by increasing the num- 
ber of rays shot. In practice, one always needs to look for 
the best compromise with computational efficiency. 

All the above ingredients have been successfully imple- 
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merited in our code CRASH which is to date one of the main 
references among RT numerical schemes used in cosmology. 
The first version of the code (Ciardi et aZ.2001, COl; in the 
following we will refer to this version as CRASHO) integrated 
the principles of ray-tracing and MC to follow the evolu- 
tion of hydrogen ionization for multiple sources under the 
assumption of a fixed temperature for the ionized gas. In 
its actual reference version (Maselli, Ferrara & Ciaxdi 2003, 
MFC03; MaselH & Ferrara 2005; in the following we will 
refer to this version as CRASHl), the code has been further 
developed by including a variety of additional physical ingre- 
dients (like helium chemistry, temperature evolution, back- 
ground radiation), and a more sophisticated algorithm for 
processes like photon absorption, reemission, as well as for 
the chemistry solver. All these new ingredients make the 
code extremely versatile in term of problems that can be 
studied. 

One of the main advantages of the CRASH implementa- 
tion consists in the straightforward way in which the mat- 
ter/radiation interactions are modeled and calculated. This 
allows to easily implement new physics in the code with- 
out the need for major changes in its structure. Illustra- 
tive examples are given by MCLya (Verhamme, Schraerer 
& Maselli 2006) and CRASHa (Pierleoni, Maselli & Ciardi, 
2007). Adopting the CRASH algorithms for ray tracing, spec- 
tra sampling, photon propagation and optical depth sam- 
pling (as described in Ciardi et a?.2001), we have imple- 
mented the physics of line resonant scattering in MCLya, a 
tool mainly dedicated to the study and modeling of Lya 
Emitters (LAEs). Moving forward in this direction, we ex- 
tended CRASH and implemented CRASHa, the single RT code 
to date which can self-consistently account for the simulta- 
neous and coupled transfer of line and ionizing continuum 
radiation. This makes of CRASHa a unique tool for study- 
ing the implications of high redshift LAE surveys on cosmic 
rcionization and 21 cm line from neutral hydrogen in the pri- 
mordial universe. More recently, we have started including 
molecular hydrogen in CRASH to study the inhomogeneous 
H2 photo-dissociation occurring at the dawn of the universe 
or in star forming regions. Once again, this requires very 
little manipulation on the main core of the code and can be 
done by complementing the existing numerical scheme with 
algorithms which model the new physics to be implemented. 

As the code is continuously being updated and opti- 
mized, in this paper we describe some major changes that 
have been introduced in the implementation described in 
MFC03, together with some test results clearly showing the 
improvements achieved in performance and accuracy. 

The rest of the paper is organized as follows. We first 
give a brief review on the basics of the CRASH implemen- 
tation in Sec. 2. In Sec. 3 we describe the updates done 
in the implementation and the new algorithms introduced, 
presenting the results from test runs designed to show the 
improvements obtained. Finally we summarize our work in 
Sec. 4. 



2 CRASH: A BRIEF DESCRIPTION 

In this section we give a very brief and general description 
of the CRASH code. A more detailed and complete report is 
given in MFC03 and also in MaselU & Ferrara (2005), where 



a refined and improved version of the implementation of the 
background radiation is described. We refer the interested 
reader to these works for more details on the implementa- 
tion. 

As discussed in the introduction, CRASH is a radiative 
transfer numerical scheme which applies to 3D Cartesian 
grids and which is based on MC ray tracing techniques that 
are used to sample the probability distribution functions 
(PDFs) of several quantities involved in the calculation, e.g. 
spectrum of the sources, emission direction, optical depth. 
The algorithm follows the propagation of the ionizing radi- 
ation through an arbitrary H/He static density field and, at 
the same time, computes (i) the evolution of the thermal 
and ionization state of the gas and (ii) the distortion of the 
ionizing radiation field due to matter-radiation interactions, 
which produce the typical RT effects of filtering, shadowing 
and self-shielding. The radiation field is described in terms of 
photons which, for computational convenience, are grouped 
into photon packets and shot through the box. Both multi- 
ple point sources, located arbitrarily in the box, and diffuse 
radiation fields {e.g. the ultraviolet background or the radi- 
ation produced by H/He recombinations) are treated. 

The energy emitted by each ionizing source is dis- 
cretized into photon packets emitted at regularly spaced 
time intervals. More specifically, the total energy radiated 
by a single source of luminosity Ls , during the total simula- 
tion time, tsim, is Es — /q°"" Ls(ts)dts. For each source. Eg 
is distributed in A^p photon packets, emitted at the source 
location. In this way, the total number of photon packets 
emitted by each source, ATp, is the main control parame- 
ter which sets both the time and the spatial resolution of a 
given run. The emission direction of each photon packet is 
assigned by MC sampling the angular PDF characteristic of 
the source. 

The propagation of the packet through the given density 
field is then followed and the impact of radiation-matter in- 
teraction on the gas properties is computed on the fly. Each 
time a photon packet pierces a cell k, the cell optical depth 
to ionizing continuum radiation, , is estimated summing 
up the contribution of the different absorbers (H i , He i , 
He II ) . As the probability for a single photon to be absorbed 
in such a coll is: 

P(r,^) = 1 - (1) 

the number of photons absorbed in the cell k is the frac- 
tion P{t!^) of packet content when it enters the cell. The 
total number of photons deposited in the cell is then dis- 
tributed among the various absorbers according to the corre- 
sponding contributions given to the total cell optical depth. 
From these quantities the discrete increments of the ioniza- 
tion fractions and temperature are calculated: Axhi, AxHei, 
AxHeii, AT and (Ar)„, where AT is the positive increment 
in temperature due to photo-heating and (AT)„ the nega- 
tive one associated to the increased number of free particles. 
The trajectory of the photon packet is followed until its pho- 
ton content is extinguished or, if periodic boundary condi- 
tions are not assumed, until it exits the simulation volume. 

The time evolution of the gas physical properties (ion- 
ization fractions and temperature) is computed solving in 
each cell the appropriate set of discretized differential equa- 
tions each time the cell is crossed by a packet. In addition to 
the impact of photo- ionization and photo-heating, that, as 
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discussed above, are included in the form of discrete contri- 
butions, the implemented chemistry network explicitly in- 
cludes the other relevant radiative processes: recombination 
and coUisional ionization for the ionization fractions and var- 
ious cooling processes for the temperature (bremsstrahlung, 
Compton cooling/heating, coUisional ionization cooling, col- 
lisional excitation cooling and recombination cooling). All 
these radiative processes are instead implemented in the 
code as continuous, namely by integrating their associated 
rate coefficients. 

We again refer the reader to the documentation men- 
tioned above for further details on the integration technique 
and on the technical implementation devised to treat mul- 
tiple sources. 



3 WHAT IS NEW IN CRASH 
3.1 Colored Packet 

In this section we describe the new implementation (in the 
following we will refer to it as CRASH2) developed to deal 
with ionizing sources with arbitrary spectra. In CRASHl, each 
point and diffuse ionizing source is assigned with an arbi- 
trary spectrum which is reproduced by MC sampling the 
given spectral energy distribution (SED) and by emitting 
a large number of monochromatic photon-packets. This im- 
plementation is highly versatile, allowing to reproduce arbi- 
trary spectral shapes in a continuum range of frequency, but 
has two main drawbacks that could affect the precision and 
efficiency of the code. First, the nature of the MC sampling 
procedure induces a poor sampling of the spectral region 
with lower intensity, which in most applications of interest is 
the hard energy tail of the spectrum. This limitation could 
be crucial in studies involving e.g. helium (He ii ) photo- 
ionization, hard radiation photo-heating (X-rays), ionization 
front (IF) profiles of hard sources like quasars, just to give 
some examples among many. Second, the monochromatic 
packet approach aggravates the well known shortcoming of 
ray-tracing techniques, which oversample the regions close 
the sources at the expense of undersampling the farthest 
ones. This introduces unphysical features in the RT results, 
such as spikes and uneven ionization fronts [e.g. see Fig. 12 
in 106). 

To overcome these problems, we have introduced a new 
algorithm in CRASH which treats the polychromatic nature 
of the ionizing radiation field by shooting "colored" pack- 
ets. Each packet is now represented by a vector of A*'^ bins, 
whose total energy Ai5 (determined as in CRASHl) is dis- 
tributed among a number of photons which populates the 
Nv frequency bins to reproduce the SED characteristic of the 
ionizing source. The number and spacing of bins can be ar- 
bitrarily chosen together with the frequency range covered. 
Typically, the spacing is chosen to be logarithmic and more 
refined at the photo- ionization thresholds, i.e. um, J^HcI and 
I'Hcii. Although we now have polychromatic photon packets 
traveling through the box, their interaction with the gas has 
not been modified significantly with respect to CRASHl. 

While propagating through the gas the colored packet 
sees different optical depths, one for each frequency bin: 

i=HI,HcI,HcII i=HI,HcI,HeII 



where the indexes i and j span over the species included 
and the frequencies sampled, respectively. The other sym- 
bols in eq.[2]have the conventional meaning, at and rii being 
the photo-ionization cross section and the number density 
of the i specie, and 51 the actual path that the packet travels 
inside the cell. Differently from the CRASHl implementation, 
here we have introduced a new algorithm which computes 
the length 51 of the ray segment within each cell (see next 
section). Once the opacity of the cell is determined, the num- 
ber of photons deposited is calculated according to the prob- 
ability for a photon to be absorbed, Pivi) = (1 — e"^'"'"'*'), 
which depends on the atomic species and the frequency. 

Summarizing, a packet is emitted from a source with 
N'^{vi) photons per frequency bin at the time tp. While trav- 
eling through the gas, photons are progressively removed 
from the packet as they ionize atoms along the path. The 
number of photons deposited in a cell is evaluated from the 
total optical depth, which gives the probability for an ab- 
sorption event: 

= ^t:;"(^.) (1 - e--'"^)) , (3) 

where c counts the number of cells crossed by the packet 
from its emission location, Aa,7 the number of photons ab- 
sorbed and At,7 the number of photons which is transmitted 
to the next cell. 

The absorption probability in the various frequency bins 
induces a modification in the spectral shape of the packet. 
Typically, the frequencies closer to the resonant photo- 
ionization thresholds are depleted first and more strongly. 
In Figure [T] we show an example of how the population of 
the frequency bins changes while the packet travels through 
a H/He gas. The differently colored curves represent the 
spectral shape of the packet in the different cells pierced, 
with the sequence from top to bottom corresponding to 
the order in which they are crossed. Note that for both 
panels each color represents the cell after a given number 
of crossings, e.g. the transitions dark blue-to-magenta and 
purple-to-green roughly correspond to the 20th and 30th 
cell crossed. The two panels show the evolution of the spec- 
tral shape of a packet emitted (i) in a completely neutral 
gas (left) and (ii) when the H ii / He ii ionization fronts are 
roughly thirty cells away from the source (right j3- It can be 
seen how the photons close to the photo-ionization thresh- 
olds are absorbed first, due to the higher cross sections, and 
how the SED of the packet becomes progressively harder as 
it moves far away from the source. It is interesting to notice 
that the spectrum evolution shown in the left panel does not 
exhibit the drop at 54.4 eV, as here the packet is emitted in 
a completely neutral gas, still devoid of He ii absorbers. 

A proper modeling of this filtering mechanism is very 
important, particularly to correctly estimate the gas tem- 
perature or the inner structure of the ionization fronts. 

Once the total number of photons deposited in a cell 
from each frequency bin Uj has been determined, NA,-y(vi), 
they are distributed among the absorber species included 
in the calculation. This is done by assigning to each ab- 



The physical configuration of the simulation analyzed here is 
the same of the CRASH/CLOUDY comparison described in Sec. 4.3, 
and the two panels correspond to packets emitted at the simula- 
tion times is = yr (left) and is = 6 X 10^ yr (right). 
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Figure 1. The curves show the evolution of a photon packet's spectrum as the paclcet travels through the simulation volume, in a 
H/He gas. From top to bottom, each curve corresponds to the spectrum at each cell in the crossing order. Le/t: The packet is emitted 
in a completely neutral gas. Right: The packet is emitted when the H II / He II ionization fronts are roughly thirty cells away from the 
source. 



sorber species a fraction of NA,-y{i^j) proportional to the 
corresponding absorption probability 

P,(/.,)(l-e-^'(^^-'), 



normalized to the total probability: 



E 

i=HI,HcI,HoII 



(4) 



(5) 



As the packet passes through the cell, it photo-ionizes a 
number of atoms of the i-th species given by: 



(6) 



with a corresponding increment in the ionization fractions 
given by: 



Nu 



(7) 



where is the singly ionized ion of the i specie (i^ £ {H ii , 
He II ,He iii } for i £{H i , He i ,He ii }), and Na the total 
number of nuclei of a given element contained in the actual 
cell. 

The cell temperature before the current packet cross- 
ing, Tc, is updated accounting for photo- heating and for the 
variation in the number of free particle which is given by 
the variation in the number density of free electrons He- 
The two effects contribute respectively with the following 
temperature increments: 

2 Y^i^^i-^i-^^J - ^^thA) 



ATn = 



2 AWe 
6 Tie 



(8) 
(9) 



where hvth,i is the ionization threshold for photoionization 
of the i-th specie, ric is the number density of free particles 



and Aric is the variation of this quantity corresponding to 
the change in the electron number density due to photo- 
ionization: 

An 

= Axe = [fh X Axmi + fhe X (Aa;Hoii + 2Aa::Hciii)] . (10) 

ric 

Each time a colored packet crosses a new cell the set of 
vectorial operations described in equation [2] through [10] is 
performed to get the key quantities: the number of photons 
which escape the c-th cell crossed from the emission loca- 



tion, Ar;^';^^(iv,) = iV, 



(c-l) 



(c) 

N\ , which accounts for radiative 



transfer, and Ax^+, AT^+ and AT„ necessary to account for 
the evolution of the gas physical state. 

The radiative processes which are not associated with 
the ionizing radiation are reproduced, as done in CRASHl, by 
integrating the corresponding temperature dependent rates 
over the elapsed time interval between the previous packet 
crossings and the current one. 



3.2 Ray Casting 

To limit the computational time, in the previous versions of 
CRASH ray casting is not implemented. Rather, the opacity 
of a single cell is evaluated using as the length traveled by 
the packet within the cell the fix value 0.56 (in units of the 
linear dimension of a cell). This number corresponds to the 
median value of the PDF for the length of randomly oriented 
paths within a cubic box of unit size (see Ciardi et aZ.2001). 

Although this approximation is acceptable for several 
applications, it does not always result in enough accuracy, 
occasionally showing spurious geometrical patterns in the 
ionization and temperature maps (see the rings evident in 
the bottom right panel of Figure [5}. For this reason, a new 
algorithm for ray casting has been implemented, which is 
similar to the Fax Voxel Trasversal Algorithm for Ray Trac- 
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ing by Amanatides & Woo (1987). Ray casting is performed 
by computing the distances that need to be traveled in or- 
der to intersect the boundaries (along the x, y and z axis) 
of the grid cell pierced by the ray. The minimum of these 
distances gives the length of the path traveled within the 
cell and the boundary intersected first determines which in- 
dex has to be incremented to get the coordinate of the next 
cell through which the ray passes. For more technical details 
about the algorithm we refer the reader to the Amanatides 
& Woo (1987) paper. 

The test cases described in Sec. 4 show the improvement 
resulting from the introduction in CRASH2 of the ray cast- 
ing calculation described above. This improvement in per- 
formance comes at the expense of computational efficiency, 
which nevertheless is not dramatically affected. The fraction 
of total computational time spent in computing ray casting 
strongly depends on the specifics of a given run, particularly 
on the physics included. In all the cases analyzed here we 
find it to be almost negligible, being at most five percent of 
the total computational time. 

3.3 Tables and Rates 

In the new version of the code we have updated the expres- 
sions adopted for the coefficient rates. The references for 
each specific coefficient can be found in Table 1 of 106. We 
have furthermore introduced the use of pre-compiled tables 
for the rate coefficients with the aim of improving efficiency. 
The temperature dependence of the various coefficient rates 
are now computed once at the beginning of a given run, 
stored in dedicated tables and read during the run. This al- 
lows to speed up the calculation significantly with respect 
to the previous implementation where the coefficients rates 
are functions evaluated at each cell crossing. Although a 
quantitative estimate of the improvement in the efficiency 
could be given, this is highly dependent on the problem at 
hand and would be significant only in specific cases. Typi- 
cally though, we find that the efficiency is increased up to 
^ 20 percent. The implementation is worked out to allow for 
adaptive tables, by assigning as initial conditions the range 
of temperatures and the refinement of sampling according 
to the physical configuration associated to a given run. This 
is done once, at the beginning of each run, with the aim of 
adapting the range and refinement of the tables to the phys- 
ical configuration of the simulation at hand: e.g. if the gas 
in the simulation has been shock heated we assign an upper 
limit on the temperature range of the order of Tmax ~ lO^K, 
which is not the case if we know that the run will process 
only photoionized gas, for which we adopt Tmax ~ lO^K. 



4 RESULTS AND PERFORMANCE 

In the following we will discuss some test cases run to study 
the performance of the new implementation described above. 

4.1 Comparison Project Test 2 

As a first test we have repeated Test 2 of the RT Com- 
parison Project (106). This test is the classical problem of 
an H II region expansion in a uniform gas around a single 
ionizing source. A steady source with a 10^ K black-body 



spectrum, emitting iV.^ = 5 x 10** ionizing photons per unit 
time (s~^), is turned on in an initially-neutral, uniform- 
density, static environment with hydrogen number density 
riH = cm~^. The computational box dimension is L 

= 6.6 kpc and the source is positioned at its corner. We 
follow the evolution of the neutral hydrogen fraction and 
of the gas temperature (initially set to 100 K), running the 
CRASHl version with lO" photon packets and the CRASH2 one 
with 2 X 10* photon packets. 

The H I fraction and temperature maps, 10 Myr after 
the source has turned on, are shown in Figure [2] for CRASHl 
(right) and CRASH2 (left) (same cuts in Figs. 11 and 12 of 
106). From a comparison between the left and right panels, 
differences are evident. More specifically, the spikes present 
in the CRASHl maps, particularly in the temperature, are al- 
most completely vanished. This is due to the colored packets 
implementation that allows for a much better sampling of 
the high energy tail of the spectrum, for the same num- 
ber of packets emitted (and thus for the same resolution). 
The weak spikes-like features which are still visible in the 
CRASH2 temperature map, are due to the noise intrinsic to 
Monte Carlo algorithms. Also, the artificial ring-like pat- 
terns clearly visible in the CRASHl maps, have completely 
disappeared thanks to the new ray casting algorithm imple- 
mented. Finally a visual inspection of the maps in Figure [2] 
reveals that the volume inside the I-fronts is both hotter and 
more highly ionized in CRASH2 than in CRASHl. 

This can be better seen in the more quantitative com- 
parison shown in Figure |3] which plots the spherically- 
averaged temperature (left panel) and ionized and neutral 
fractions (right panel) as a function of the distance from 
the source 100 Myr after the source has been turned on. 
CRASH2 (solid red lines) produces a higher temperature and 
a lower neutral fraction compared to CRASHl (dashed blue 
lines). This result, which brings a better agreement between 
CRASH and other codes {e.g. C^-ray, whose results are shown 
for comparison by the dashed-dotted green curves), is due 
once again the improved modeling of the filtering which is 
achieved thanks to a better sampling of the high energy tail 
of the spectrum: photons in this region of the SED can heat 
more and ionize less the gas with respect to lower energy 
photons which are oversampled in CRASHl. Despite this, as 
a result of the increased temperature which inhibits recom- 
binations, the ionization fraction within the ionized sphere 
results higher in CRASH2. One can also see that the ioniza- 
tion front lies slightly closer to the ionizing source in the 
new implementation: this happens because, once the lumi- 
nosity is assigned as energy emitted in the unit time, the 
oversampling of lower energy part of the SED (which afi'ects 
the CRASHl implementation) results in an higher number 
of ionizing photons. Also, we still obtain a wider ionization 
front with respect to other codes (see Fig. 17 in 106), which 
results from the better sampling of the spectrum in the fre- 
quency space. 

4.2 Comparison Project Test 4 

In this section we describe the results from the Test 4 of 
the RT Comparison Project. In this case, the propagation 
of ionization fronts from multiple sources in a static cosmo- 
logical density field is followed. The initial conditions are 
provided by a time-slice (at redshift z = 9) from a cosmo- 
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Figure 2. Test 2 of the RT Comparison Project (H II region expansion in a uniform H gas with varying temperature). The upper (lower) 
panels show maps of the H I fraction (temperature) cut through the simulation volume at coordinate z = and time t = 10 Myr. Left: 
CRASH2; Right: CRASHl. 



logical N-body and gas-dynamic simulation. The simulation 
box size is Q.^h~^ comoving Mpc. The gas is assumed to be 
initially neutral and the temperature is initialized at 100 K 
everywhere in the box. The ionizing sources are chosen so 
as to correspond to the 16 most massive halos in the box, 
which are identified with a friends-of-friends algorithm. We 
assume that the sources have a black-body spectrum with 
effective temperature Te// = 10^ K. The ionizing photon 
production rate for each source is constant and is assigned 
assuming that each source lives for 3 Myr and emits 250 
ionizing photons per atom during its lifetime. For simplicity, 
all sources are assumed to switch on at the same time. The 
boundary conditions are transmissive {i.e. photons leaving 
the computational box are lost, rather than coming back in 
as for periodic boundary conditions). The maps in Figure |4] 
show the H i fraction (upper panels) and temperature (lower 
panels) cuts through the simulation volume ai z — Zbox/'2 
and at time 0.05 Myr (same cuts in Fig. 31 and Fig. 32 of 
106). Both CRASHl (right panels) and CRASH2 (left panels) 
outputs are shown in the figure. The differences between 
the two implementations are strikingly evident: as already 
found in the Test 2 above, the ionization fronts are much 
larger in the CRASH2 implementation which accounts prop- 



erly for the higher energy photons associated with larger 
mean free paths. The better sampling of the harder tail part 
of the spectrum results also in higher temperature and lower 
neutral hydrogen fractions as already discussed above. The 
new ray casting implementation corrects furthermore for the 
artificial patterns visible close-by the ionization fronts in the 
CRASHl. The difference in the internal structure of the ion- 
ized region can be better appreciated in Figure [S] where the 
histograms of the temperature at time t = 0.05 Myr ob- 
tained with the two implementations of the code are com- 
pared: solid (dotted) hne is for CRASH2 (CRASHl). With the 
colored packets implementation we find on average higher 
temperatures. The fact that regions at temperature above 
10^ K are more densely populated reflects the higher tem- 
perature of ionized regions produced by the contribution of 
harder photon packets which are undersampled in CRASHl. 
At the same time the lack of the peak found at low tem- 
peratures in CRASHl is due to the wider ionization fronts 
obtained with CRASH2, which extend much further from the 
sources and affect almost the whole computational volume. 
As in the previous test, with CRASH2 we found a better 
agreement with C^-ray results shown in the figure by the 
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Figure 3. Test 2 of the RT Comparison Project (H II region expansion in an uniform H gas with varying temperature). The left (right) 
panel shows the spherically-averaged temperature (ionized and neutral H fraction) as a function of the distance from the source, 100 Myr 
after the source has turned on. Solid line: CRASH2; Dashed line: CRASHl; Dashed-Dotted line: C^-RAY. 
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4.3 CRASH/CLOUDY Comparison 

Finally, similarly to what was done in MFC03, we compare 
our code, with the publicly available ID RT code CL0UDY9^. 
With this test we seek to check for the performance of the 
new algorithm in dealing with helium physics, which has not 
been included in the comparison project tests. We consider 
the case of a point source with a black body spectrum at 
T = 6 X 10* K and luminosity L = 10^* erg s~^, ionizing a 
uniform medium with density n — 1 cm~^ in a gas composed 
by hydrogen (90 percent in number density) and helium. Fig- 
ure|6]shows the profiles at equilibrium for temperature, neu- 
tral and ionization fractions. Compared to the earlier results 
of CRASHl (see Fig. 4 in MFC03), the new version of the code 
provides an overall much better agreement with CL0UDY94. 
In particular the neutral hydrogen fraction inner profile fits 
extremely well and with reduced dispersion the CL0UDY94 
predictions. The same is true for the neutral helium frac- 
tion profile which is however slightly higher in CRASH2 than 
in CLDUDY94. Nevertheless, both the temperature and the 
He 11 profiles are slightly lower than the ones derived with 
CL0UDY94 and the ionization front extends somewhat further 
in CL0UDY94. The two behaviors are likely to be correlated 
and might be due to differences in the expressions adopted 
for the coefficient rates associated to helium together with 
the fact that CRASH neglects the effects of heat conduction. 



C 




Figure 5. Test 4 of the RT Comparison Project (reionization 
of a cosmological density field): histogram for the temperature 
at time t = 0.05 Myr. Solid line: CRASH2; Dashed line: CRASHl; 
Dashed-Dotted line: C^-RAY. 



It should be noted that in Fig. 37 of 106 the purple and green 
lines have been erroneously exchanged. 
^ http:/ /nimbus. pa. uky.edu/cloudy 



5 SUMMARY 

In this paper we introduce the new version of the radia- 
tive transfer code CRASH. The details of the code imple- 
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Figure 4. Test 4 of the RT Comparison Project (reionization of a cosmological density field). The upper (lower) panels show maps of 
the H I fraction (temperature) cut through the simulation volume at coordinate z = zi,ax/2 and time t = 0.05 Myr. Left: CRASH2; Right: 
CRASHl. 



mentation have been previously described in MFC03. Here 
we present the new elements introduced to improve per- 
formance and efficiency. The main new ingredients imple- 
mented in CRASH2 are: 

• the introduction of colored photon packets, which 
brings a significant improvement in the sampling of source 
spectra both in space and frequency; 

• a new model for propagating the photon packets 
through the Cartesian grid, which includes a proper ray cast- 
ing algorithm; 

• the use of pre-compiled adaptive tables for the coeffi- 
cient rates which result in more efficient computing. 

The CRASH2 implementation has been tested and compared 
with the earlier versions of the code by performing a set of 
reference tests taken from the RT codes comparison project 
(106) and from MFC03. CRASH2 is more accurate both in 
the calculation of the temperature and the ionization frac- 
tions, thanks to the improved sampling of the hard tail part 
of ionizing sources spectra, which is undersampled in the 
monochromatic packet implementations. Furthermore, the 
better sampling of hard radiation together with the new ray 
casting algorithm corrects for spurious geometrical features 



like spikes in Stromgren sphere boundaries or artificial geo- 
metrical pattern that are found in the CRASHl results. 
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Figure 6. Comparison between tlie equilibrium configurations obtained by CRASH2 (points) and CLQUDY94 (solid lines) for different 
physical quantities, as a function of distance from the point source in cell units. The errors refer to the CL0UDY94 results and are evaluated 
as explained in the text. 
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